Your friend has been patiently waiting for a design for an automatic dam control system. They've provided you specifications for the design and now you are going to have to use your knowledge of transfer functions, including steady state and transient behavior, to design a closed-loop controller to regulate the dam's turbine speed.
In your last notebook you were introduced to the "canonical" form of a system under feedback control, which is shown in the diagram below:
In the diagram above, as labeled, the system has the following key pieces:
We saw in the last reading that the canonical feedback control system has a closed-loop transfer function:
\begin{equation} \frac{y(s)}{u(s)} = \frac{CG}{1+CGH} \end{equation}But what does this really mean? Note that the controller transfer function, chosen by the engineer, shows up in both the numerator and the denominator of the closed loop transfer function. This means that the addition of the controller and feedback to the plant changes the dynamics of the plant output. This is precisely what we want, if our goal is to improve system performance. In simple terms, $C(s)$ is a transfer function that "makes decisions" about how hard to "push" the plant based on the difference between where the system is at any given instant ($y(s)H(s)$), and where the operator wants the system to be ($r(s)$).
In this course, we will focus on variants of the PID or "proportional-integral-derivative" controller, which has the following general transfer function:
\begin{equation} C(s) = \frac{u(s)}{e(s)}= K_p + K_i\frac{1}{s} + K_ds \end{equation}The PID controller responds to the error proportionally (P-control), and/or the error's derivative proportionally (D-control), and/or the error's integral proportionally (I-control). PID control almost always includes a P-control element in the form of $K_p$, but often, $K_i$ and/or $K_d$ can be set to zero depending on the needs of the engineer and/or the plant. This week, we will focus on how variants of this general design work for simple first and second order systems. We'll use simplified examples and the so-called "direct method" of determining closed-loop system performance (see, for more information, chapters 3 and 4 of Schaum's Outline of Feedback and Control Systems or Sections 13.1-13.3 of Dynamic Modeling and Control of Engineering Systems) to build an understanding of how each type of controller works before diving into some more involved design tools that will allow us to apply these controllers effectively to more complex systems.
When synthesizing P(I)(D) controllers, we can have a number of different design goals. Some of these are listed below:
Of course, following this design process, you should try to validate your design by implementing it on the real system and see if it performs as expected.
Much of what we'll discuss roughly follows the material in Chapter 10 of Schaum's Outline of Feedback and Control Systems. We will deal with how to handle most of these design goals as needed, but it is helpful to know where we are going. To begin our journey, let's take a closer look at so-called "P-control," which is a special case of PID control with $K_d$ and $K_i$ set to 0.
Proportional control is the simplest form of the PID feedback control paradigm. This controller 'pushes' on the plant proportionally to the error between the reference input and the plant output. As the system gets closer to the desired value of the output, the controller pushes less and less hard. You might think of this controller like a "spring" that drives the plant to the desired output. For a proportional controller, the controller transfer function is a constant $K_p$. Formally, we write:
\begin{equation} \frac{u(s)}{e(s)}=C(s) = K_p \end{equation}Knowing the units of your controller's gain(s) is extremely important. Often $e(s)$ and $u(s)$ will often have different units. This means that the value of $K_p$ will assume the units needed to make a dimensional analysis work out.
For example, if $e(s)$ is a distance error measured in meters, and $u(s)$ is a voltage input to a motor measured in Volts, then the controller's proportional constant $k_p$ must have units of $\frac{V}{m}$.
Fixed Assumptions
Consider a plant $P(s) = \frac{1}{s+1}$. Design a proportional controller to make this system track a step input, and achieve steady-state in less than 1 second. The system uses a sensor that measures the plant's output $y(s)$ perfectly with a gain of $1$. After designing your controller, find the maximum controller output $u(s)$ needed for a unit step input, and find the system's steady-state error between reference input and system output.
Process
To decipher the design goal, note that we often assume that a first-order system reaches steady state in roughly 5 time constants. Therefore, we can say that we'd like our system to have a time constant of 0.2 seconds or less.
Because the problem also states that the sensing in the system is 'perfect,' we can say that $H(s)=1$ as this means the output will feed back to the system without any modifications.
We will chose a proportional controller to this system and use the canonical form of the closed loop transfer function developed from block diagram simplification. Substituting the proprotional control law and our plant into this form of $\frac{y(s)}{r(s)}$ yields
\begin{equation} \frac{y(s)}{r(s)} = \frac{K_p\frac{1}{s+1}}{1+K_p\frac{1}{s+1}} \end{equation}Clearing the fractions in the numerator and denominator by multiplying the transfer function by $\frac{s+1}{s+1}$ yields the final equation for the system's closed loop transfer function:
\begin{equation} G_{cl}(s) = \frac{y(s)}{r(s)} = \frac{K_p}{s+1+K_p} \end{equation}Now, we have to determine what value of $K_p$ will allow us to reach our design goal of $\tau<0.2s$.
Note that a first order system response can be represented generically as $\frac{A}{s+a}$. Examining the denominator, we can use it to write the characteristic equation of the system and determine the root of the equation which is the system eigenvalue.
\begin{equation} s+a=0, s=-a \end{equation}The time constant can be determined from $\tau=-\frac{1}{a}$, where $a$ is the system eigenvalue.
In the case of our controller, we can meet the design constraints by having a time constant of 0.2 seconds so we want a system with $a=-5$. Solving the system's closed-loop characteristic equation means solving:
\begin{equation} s+1+K_p=0 \end{equation}This yields a system eigenvalue $a=-1-K_p$, so for our example, we know that $K_p=4$ will satisfy our design requirement. Let's simulate the closed-loop system's response to a step input on $r(s)$.
s = tf('s');
t = 0:.001:1.5;%go a little past where we said the system should be at steady state
G1 = 1/(s+1);
Kp = 4;
Gcl = G1*Kp/(1+G1*Kp);
[ysim,tsim]=step(Gcl,t);%simulate without plotting so we store variables.
figure
plot(tsim,ysim,'k')
xlabel('Time (s)')
ylabel('System Output (?)')
As you can see, our controlled system does a pretty good job of getting to steady-state by 1 second, but there's a problem... one of our design goals was for the system to track the reference input $r(s)$. Because $r(s)$ was a unit step, we might have expected that our system would reach a steady-state value of 1, but it doesn't! To show how the system's error progresses over time, we can simply take the constant value of $r(s)$ (which in this case was just 1) and subtract the system's output under closed loop control from $r(s)$, since the error signal $e(s)=r(s)-y(s)H(s)$. Alternatively, we could use some block diagram algebra to find the transfer function $\frac{e(s)}{r(s)}$ directly. I'll leave the details of this to you to try (highly recommended for practice!), but if we do the reduction, we find that the transfer function $\frac{e(s)}{r(s)}$ for our controlled system is:
\begin{equation} \frac{e(s)}{r(s)} = \frac{s+1}{s+5} \end{equation}Applying the final value theorem can get us an answer about how much steady-state error we will have in response to a unit step input:
\begin{equation} \lim_{t\rightarrow \infty} e(t) = \lim_{s\rightarrow 0} s e(s) = \lim_{s\rightarrow 0} s\frac{1}{s}\frac{s+1}{s+5} = 0.2 \end{equation}So it appears that our controller, while it does meet the requirement to reach steady state within 1 second, doesn't do a great job of tracking our reference input! Unless we're ok with $20\%$ tracking error even for a simple step input, we'll need something else to do better. Let's plot the system's error response using MATLAB/Octave:
%create our TF from e to y
G_e_y = (s+1)/(s+5);
%now see how it responds to a step input
[yey,tey] = step(G_e_y,t);
figure
plot(t,yey,'k')
xlabel('Time (s)')
ylabel('tracking error')
It appears that the error acts as we'd expect. The error signal could also have been found by simply subtracting $y$ from $r$ given our original step response. To see how much system input would have been required to achieve this step response from our controller, we need to plot the step response of the transfer function $\frac{u(s)}{r(s)}$, or in this case, we can simply multiply $e(s)$ by $K_p$. This implies that the maximum controller output needed for the unit step is $4$ "units." These units could be voltage if our controller drove an electric motor, pump, solenoid, or similar electromechanical device, for instance.
In this exercise, you will use what you now know about block diagrams, and combine it with your linearized model relating penstock valve position to turbine speed to design a proportional control system to regulate turbine speed by modulating the valve position.
This means that our design will only be valid for small requests in changing turbine speed from the turbie's NOP.
A simulator is provided for you below.
%%html
<iframe id="inlineFrameExample" title="Inline Frame Example" width="1000" height="900" src="https://workbench.lafayette.edu/~brownaa/ME480/damsim/damsim_w5m.html"> </iframe>
To confirm that your model fits the turbine's behavior near an NOP, re-produce the validation step from Week 4 Monday. If your model did not work well then, this is your chance to improve it. If (and only if) you do make model changes, paste them in the markdown cell below. If you do not make changes, simply reproduce your results from Week 4 Monday, section 2.2.2 in the code cell below. You must show model agreement with data on one set of axes for the following experiment:
Using the simulator, set the initial condition by beginning with the dam overflowing and the valve open at 40% before data collection begins. Allow the turbine to come to a steady state speed. Once the system is at steady state, the user presses a button, the valve percentage should be increased to 50% while data are collected for 5 seconds. Be sure the electrical system is disconnected for this test
no changes were made to my linearized model from Week 4 Monday
data = load('damdata_lin.txt');
Pin = data(1,3)*1000;
tnl = data(:,1);
ynl = data(:,4)*2*pi*(1/60);
%Set up parameters of the models
J = 4250;
b = 2000;
D = 1;
Rnom = 1.936E4;
%simulation parameters
simtime = 5;
dt = .001;
%simulation time vector-- we can make it common between both nonlinear and linear models
t = 0:dt:simtime;
%set up input situation
ubar = 40;
du = 10;
%set up the input vector for TOTAL input u = ubar + du
duvec = du*ones(size(t));
u = (ubar+du)*ones(size(t));
%set up the initial conditions. Use the NOP equations to set them!
ybar = (D*Pin(1))/(b+(D^2*Rnom)-(0.01*D^2*Rnom*ubar));
%set up variables for the linear simulation. INCREMENTAL variables!!
Dyl = zeros(size(t)); % y for linear
Dydl = zeros(size(t)); %ydot for linear
%set up the initial conditions. Use the NOP equations to set them!
% For linearized, we are using Delta y, so the IC is 0!!
Dyl(1) = 0;
Dydl(1) = 0;
%perform Euler Integration for the LINEARIZED model
for k = 2:length(t)
%highest derivative is a function of the others
Dydl(k) = (1/J)*((-b-(D^2*Rnom)+(0.01*D^2*Rnom*ubar))*Dyl(k-1)+(0.01*D^2*Rnom*ybar*duvec(k-1)));
%y can be integrated using the previous value of ydot.
Dyl(k) = Dyl(k-1) + dt*Dydl(k);
end
%shift linearized model up to start at NOP!
yl = Dyl + ybar;
% plot the results to see agreement!!
hold on;
plot(tnl,ynl,'.')
plot(t,yl)
xlabel('Time (s)')
ylabel('Turbine Angular Speed (rad/s)')
legend('Nonlinear Data','Linearized Model','location','southeast')
hold off;
Design a P-controller for the dam system. You know that the "correct" valve position is somewhere around 40% to achieve the total turbine speed $\Omega$ you want the dam to run at. However, you want to be more precise than this, so you implement a proportional feedback controller to help regulate the turbine's speed.
Your P-controller should modulate the change in valve position above or below 40%, which we called $\Delta u_v$ in Week 4 Monday, to slightly open or close the valve to achieve a desired turbine speed.
Your controller design goal should be to achieve a steady-state error between desired and true change in turbine speed, or $e=\Delta \Omega_d - \Delta \Omega$ of less than 20% of the requested change for a step change in desired turbine speed of 2 rad/s. All changes should occur relative to an NOP defined by a nominal valve position of $\bar{u}_v = 40\%$.
Present your hand-calculations to achieve this controller design in the markdown cell below.
In addition to the variables such as VALVE, VALVEPERCENT, and REC that you have become accustomed to using in the simulator, you now also have the ability to access the turbine's speed in radians/second by calling the variable OMEGARS. You also now have access to a global variable called Kp, which is initially set to 0. You can set Kp to any value you like inside your simulator code. Finally, you now have access to a global variable called OMEGADES, which is initially set to 0. You can set this value to anything you like to represent the desired speed for your turbine in rad/s.
These new variables will allow you to implement closed-loop control by setting VALVEPERCENT according to the current turbine speed!
The data file's structure is the same as it was for your last assignment, except that VALVEPERCENT has been added as a 7th column.
| column | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|
| data | time | gauge height (ft) | Pressure at penstock valve (kPa) | turbine speed (RPM) | generator volts | house volts | valve percent |
In the markdown cell below, update your data collection code to implement proportional control on valve position according to the current turbine speed. The operation should be as follows:
OMEGADES and the actual turbine speed OMEGARS.You should only need to update your block 4 code, so no new FSM design is required-- just a direct modification to the code will be sufficient. Properly format and paste your updated code below.
//BP represents pressing X1
//elapsed1 represents the current time
//lastStoredTime represents the last time the X1 was pressed and recording began
//V1 represents the waiting state
//V2 represents the recording state
//V3 represents the transition from waiting to recording
//V4 represents the transition from recording to waiting
//V5 represents the latch on the waiting state
//V6 reptresents the latch on the recording state
//VALVEPERCENT represents the percent that the valve is opened
//Y1 represents the light labeled Y1
//REC represents the recording of data
//OMEGADES represents the desired turbine angular speed (in this case it is 2 rad/s above the steady state at 40% valve opening)
//OMEGARS represents the actual angular speed of the turbine
//Kp represents the P-controller transfer function, which is just a constant calculated above
BP = X1;
elapsed1 = millis();
Kp = 36.2738;
if (V3) {
lastStoredTime1 = millis();
}
//BLOCK 2
V3 = V1&&BP;
V4 = V2&&(elapsed1-lastStoredTime1>=5000);
V5 = V1&&!BP;
V6 = V2&&(elapsed1-lastStoredTime1<5000);
//BLOCK 3
V1 = V4||V5||SP0;
V2 = V3||V6;
//BLOCK 4
if (V2 == 1){ // this changed from last assignment so that the valve percentage can be changed
OMEGADES = (10.987+2);
VALVEPERCENT = 40+(Kp*(OMEGADES-OMEGARS));
} else {
VALVEPERCENT = 40;
}
Y1 = V2;
VALVE = 1; // this is changed from last assignment so that the valve begins open
REC = V2;
Using the updated control code above, implement your control law on the dam. Collect data with the valve position initially open to $\bar{u}_v = 40\%$. When X1 is pressed, the desired speed should be increased by 2.0 rad/s over the turbine steady state angular velocity at the 40% valve setting.
Then, simulate your closed-loop system's behavior in response to the step change in desired turbine speed. You can use have MATLAB/Octave simulate the responses using your transfer functions and the step() command. Plot the simulated closed-loop response and your collected data on one set of axes.
data = load('damdata_pcontroller2.txt'); % loading collected data
t = data(:,1);
y = data(:,4)*2*pi*(1/60);
%Set up parameters of the model
J = 4250;
b = 2000;
D = 1;
Rnom = 1.936E4;
omegabar = 10.987; % steady state angular speed (from last week's assignment)
ubar = 40;
% P Controller Model...
s = tf('s');
time = 0:.001:5;%go a little past where we said the system should be at steady state
G1 = (0.01*D^2*Rnom*omegabar)/(J*s+(b+(D^2*Rnom)-(0.01*D^2*Rnom*ubar)));
Kp =36.2738;
Gcl = G1*Kp/(1+G1*Kp);
[ysim,tsim]=step(Gcl,time);
ysim= (2*ysim)+omegabar; % multiplying by 2 because we're increasing by 2 rad/s and adding turbine's previous steady state value when valve was open 40%
% plotting...
figure
hold on;
plot(tsim,ysim,'k')
plot(t,y,'r.')
plot(t,ones(length(t))*12.987,'k--')
xlabel('Time (s)')
ylabel('Turbine Angular Speed (rad/s)')
title('Comparison of P-Controller Model to Experimental Data from Dam System')
legend('P-Controller Model','Experimental','Desired','location','southeast')
hold off;
In the markdown cells below, evaluate your model. Answer the following two questions.
The controller does satisfy the design requirement. Using the method described above, a value for Kp was determined such that the steady state error was below 20% of the desired turbine speed (in this case the steady state error was chosen to be 15%). The desired turbine speed in this case was 2 rad/s above the steady state turbine speed when the valve is open 40%. The desired turbine speed for this case is 12.987 rad/s and is represented by the dashed black line on the graph above. To be within the 20% steady state error region, the final steady state turbine speed produced by both experimental data and simulated model should be at least 12.587 (12.987 rad/s - (0.2*2 rad/s)). The final steady state turbine speed achieved by the P-controller model and the experimental data was approximately 12.7 rad/s, which is well within the design requirement of less than 20% steady state error.
The P-controller model predicts the actual behavior of the dam quite well, as they reach the same steady state value. The only deviation is that the actual behavior seems to have a faster rise time than what the P-controller model predicts. The actual data also shows that the steady state turbine speed is initially overpredicted and later the valve percent will decrease in order to reach the proper steady state turbine value. This is most likely due to the fact that when given the initial step at t=0+, the error (omega desired - omega actual) is very large, and thus the amount that the valve is opened is intially very large. As this error decreases and the actual turbine speed approaches the desired, the change in valve percent is much less drastic. I think there is almost somewhat of a delay as within the feedback loop as the system is constantly having to adjust the percent the valve is open based on the difference between the actual and desired turbine speeds.